function D=osl_forward_model(S)

% D=osl_forward_model(S)
% 
% Required inputs:
% S.D: SPM MEG object filename
% S.mri: structural MRI nii file name (set S.mri=[] or '' to use template
% structural)
% S.useheadshape: set to 0 or 1 to indicated if the headshape points should
% be used in the registration
%
% Optional:
% S.forward_meg: e.g. set to 'Single Shell' or 'MEG Local Spheres' (default)
% S.fid_label.nasion: default: 'Nasion'
% S.fid_label.lpa: default: 'LPA';
% S.fid_label.rpa: default: 'RPA';
%
% Check the output using:
% spm_eeg_inv_checkmeshes(D);
% spm_eeg_inv_checkdatareg(D);
% spm_eeg_inv_checkforward(D, val);
%
% MWW 2012

global OSLDIR;

try, useheadshape=S.useheadshape; catch, error('S.useheadshape not specified'); end;
try, spmfilename=S.D; catch, error('S.D not specified'); end;
try, mri=S.mri; catch, error('S.mri not specified'); end;
try, correct_grads=S.correct_grads; catch, correct_grads=1; end;

try, forward_meg=S.forward_meg; catch, disp('Using Single Shell MEG head model'); forward_meg='Single Shell'; end;
try, fid_label=S.fid_label; catch, warning('Using default Neuromag fid labels'); fid_label.nasion='Nasion'; fid_label.lpa='LPA'; fid_label.rpa='RPA'; end; % To see what these should be look at: D.fiducials; fidnew.fid.label    

if(strcmp(mri,''))
    mri=[OSLDIR '/std_masks/MNI152_T1_1mm.nii'];
    
    warning(['Using ' mri ' as structural']);
end;

spmfilename_orig=spmfilename;
D = spm_eeg_load(spmfilename);

copy_fm_flag = 0;

if isfield(D,'raw_tra_matrix') && strcmp(forward_meg, 'MEG Local Spheres');

    S_copy         = [];
    S_copy.D       = D;
    S_copy.newname = fullfile(D.path, ['temp_' D.fname]);
    Dnew=spm_eeg_copy(S_copy);
   
    grad=D.sensors('MEG');
    
    %%%%%%%%%%%%%%%
    grad.tra=D.raw_tra_matrix;
    Dnew = sensors(Dnew,'MEG',grad);
       
    %gradraw=D.raw_sensors;
    %gradraw.label=grad.label;    
    %Dnew = sensors(Dnew,'MEG',gradraw);    
    %%%%%%%%%%%%%%%%
    
    Dnew.save;
    
    spmfilename=fullfile(Dnew.path,Dnew.fname);
    copy_fm_flag=1;   
    
elseif ~isfield(D,'raw_tra_matrix') && strcmp(forward_meg, 'MEG Local Spheres');
    error('Raw tra matrix not detected. MEG Local Spheres forward model may default to Single Sphere forward model unpredictably.');
end

matlabbatch{1}.spm.meeg.source.headmodel.D = {spmfilename};
matlabbatch{1}.spm.meeg.source.headmodel.val = 1;
matlabbatch{1}.spm.meeg.source.headmodel.comment = '';
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshes.mri = {[mri ',1']};
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshres = 2;
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(1).fidname = fid_label.nasion;
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(1).specification.select = 'nas';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(2).fidname = fid_label.lpa;
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(2).specification.select = 'lpa';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(3).fidname = fid_label.rpa;
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(3).specification.select = 'rpa';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.useheadshape = 1;
matlabbatch{1}.spm.meeg.source.headmodel.forward.eeg = 'EEG BEM';
matlabbatch{1}.spm.meeg.source.headmodel.forward.meg = forward_meg;

spm_jobman('run', matlabbatch);

chant=D.chantype;

if copy_fm_flag
    Dnew = spm_eeg_load(spmfilename);
    D.inv{1}=Dnew.inv{1};
    D.save;
    Dnew.delete;
else
    D = spm_eeg_load(spmfilename);
end


if sum(ismember(unique(chant),'MEGGRAD')), 
    
    disp('Using CTF data.');
    
elseif sum(ismember(unique(chant),'MEGPLANAR')) && ~sum(ismember(unique(chant),'MEGPCACOMP'))

    disp('Using Elekta Neuromag 306 data.');

    if isfield(S,'neuromag_planar_baseline_correction') && ~strcmp(S.neuromag_planar_baseline_correction,'none'),        
        D=osl_neuromag_grad_baseline_correction(spmfilename_orig,S.neuromag_planar_baseline_correction);
    end;
    
elseif sum(ismember(unique(chant),'EEG')), % added by DM
    
    disp('Using EEG data.');
    
end
